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Ill  ABSTRACT  1 

This  report  describes  FKCOMB,  a  fast  general-purpose 
array  processing  program  which  employs  f requency-wavenumber 
analysis.  A  quasi-online  long-period  array  processing 
package  has  been  designed  around  FKCOMB  for  use  at  the 
Seismic  Array  Analysis  Center.  The  input  parameters  required 
for  operating  the  program  are  described. 
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INTRODUCTION 


Frequency-wavenumber  (f-k)  spectral  estimation  is 
a  powerful  technique  for  signal  detection  and  waveform 
analysis  of  digitally  recorded  array  data.  The  f-k 
spectrum  of  a  given  segment  of  array  output  is  the 
squared  modulus  of  the  multidimensional  Fourier  transform 
of  the  data  in  time  and  space.  The  f-k  spatial  repre¬ 
sentation  of  a  propagating  wave  is  shown  schematically 
in  Figure  1,  Using  discrete  Fourier  analysis  in  the 
time  dimension,  the  representation  can  be  thought  of  as 
a  series  of  layers  normal  to  the  frequency  axis,  each 
layer  representing  the  wavenumber  plane  at  a  given 
frequency.  The  wave  is  thus  represented  as  a  series  of 
power  maxima  in  the  layers,  and  the  locus  of  these  maxima 
is  determined  by  the  phase  velocity  of  the  wave  (for 
discussion,  Smart  1971). 

The  advantage  of  this  process  is  that  propagating 
wave  components  are  easily  recognized  and  separated 
from  one  another,  subject  to  the  limitations  imposed  by 
the  array  geometry,  sensor  weighting,  and  the  type  of 
spectral  smoothing  employed.  In  essence,  f-k  analysis 
is  beamforming  in  the  frequency  domain.  The  method 
takes  advantage  of  the  fact  that  the  signal-to-noise 
ratio  varies  with  frequency ,  so  the  beamforming 
is  done  frequency  by  frequency.  Also,  by  staying 
in  the  frequency  domain  a  great  many  beams  can  be 
examined  rapidly,  the  number  being  limited  only  by  the 
resolution  cell  of  the  array  response.  In  practice  this 
means  that  the  azimuth  and  velocity  of  a  signal  need 
not  be  assumed:  one  merely  accepts  the  beam  with 
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maximum  power.  This  fact  is  important  for  signals  such 
as  long-period  seismic  surface  waves,  which  not  only 
are  dispersive,  i.e., their  phase  velocity  varies  with 
frequency  but  whose  arrival  azimuth  may  also  vary  with 
frequency  due  to  lateral  inhomogeneities  in  the  earth. 

FKCOMB  is  a  fast  f-k  analysis  program  which  was 
used  in  an  automatic  processing  system  for  microbaro¬ 
graph  array  data  (Smart  and  Flinn,  1971).  It  has  been 
incorporated  into  the  SAAC  long-period  processing 
package  for  quasi-online  processing  of  ALPA,  LASA,  and 
NORSAR  long-period  seismic  data  (Mack,  1972). 


GENERAL  DESCRIPTION  OF  FKCOMB 


FKCOMB  is  a  software  system  which  uses  f-k  analysis 
for  continuous  processing  of  time-varying  data  from 
arrays  of  sensors.  The  output  is  in  the  form  of  a 
bulletin  which  lists  signal  detections  and  gives  their 
phase  velocity,  arrival  azimuth,  signal  power,  signal- 
to-noise  ratio,  and  F  statistic  (related  to  signal-to- 
noise  ratio;  see  Smart  and  Flinn,  1971)  as  a  function 
of  frequency  and  arrival  time.  The  processor  is  shown 
schematically  in  Figure  2.  An  initial  time  window  is 
specified.  The  multichannel  data  in  this  window  are 
deglitched  and  automatically  edited  to  delete  dead  or 
noisy  traces.  The  data  are  then  Fourier  transformed  in 
time  and  space.  Maxima  of  power  in  three-dimensional 
f-k  space  are  automatically  picked,  and  if  these  maxima 
exceed  a  specified  signal- to-noise  detection  threshold, 
and  are  within  a  specified  phase  velocity  range,  the 
maxima  are  listed  together  with  the  corresponding 
arrival  azimuth,  phase  velocity,  period,  signal  power 
estimate,  signal- to-noise  ratio,  and  F  statistic.  Two- 
dimensional  maxima  also  occur;  these  are  places  which 
are  maximum  within  a  given  wavenumber  plane  but  not 
along  the  frequency  axis.  If  such  two-dimensional 
maxima  satisfy  the  specified  signal- to-noise  ratio  and 
phase  velocity  criteria,  and  if  the  corresponding 
approximation  to  group  velocity  (see  below)  is  reason¬ 
able,  then  the  maxima  are  also  listed  by  the  processor. 

An  option  is  available  at  this  point  to  detect 
smaller  signals  which  might  be  masked  by  the  detected 
signals.  In  this  option,  the  maximum  power  peaks  and 
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their  associated  sidelobes  are  removed  by  a  filtering 
process  referred  to  as  'stripping',  and  the  secondary 
maxima  are  picked  from  the  stripped  data.  If  any  of 
these  secondary  maxima  satisfy  the  original  detection 
criteria,  their  azimuth,  etc.,  are  also  listed. 

The  time  window  is  then  shifted  forward  or  back¬ 
ward  in  time  by  any  specified  amount,  and  the  entire 
process  (including  editing)  is  repeated.  A  sample  of 
the  output  listing  is  shown  in  Figure  3. 

Data  editing 

After  the  data  are  despiked,  the  variance  of  each 
channel  in  the  selected  time  window  is  calculated  and 
the  mean  variance  across  channels  is  computed.  Dead 
channels  have  very  small  variance  after  the  mean  is 
removed,  and  excessively  noisy  channels  have  a  large 
variance.  If  a  given  channel  variance  is  greater  or 
smaller  than  the  mean  by  a  given  factor  (input  para¬ 
meter  SLOP),  that  channel  is  discarded,  the  mean  of 
the  remaining  channels  is  computed,  and  this  process 
is  repeated  until  the  variances  of  the  remaining 
channels  all  lie  within  the  specified  bounds.  The 
discarded  channel  numbers  are  listed  in  the  output 
bulletin  for  each  time  window. 

The  variances  of  the  broad-band  time  domain  data 
are  used  for  editing  rather  than  the  variances  of  the 
individual  frequency  components  for  the  following 
reason:  two  signals  interfering  in  the  area  of  the 
array  yield  spectral  nulls  in  the  data  from  certain 
sensors  at  certain  frequencies.  Hence  a  large  variation 
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in  the  variance  of  a  particular  frequency  component  can 
be  due  to  a  physically  reasonable  situation  and  not 
necessarily  due  to  noisy  or  dead  traces „  Variations 
caused  by  wave  interference  are  greatly  reduced  when 
the  broad-band  data  are  used  instead. 


Time  gaps  and  time  reversals 

The  multichannel  data  are  normally  analyzed  as  a 
series  of  overlapping  time  windows,  the  length  of  the 
windows  and  the  amount  of  overlap  being  input  parameters 
to  the  program.  If  a  gap  of  one  or  two  data  samples 
occurs  within  a  particular  time  window,  then  the  last 
good  value  is  repeated  to  fill  up  the  gap,  and  the 
analysis  continues.  If  the  gap  is  more  than  two  samples 
long,  then  the  processor  analyzes  the  last  full  window 
up  to  the  data  gap  and  the  first  full  window  after  the 
gap,  continuing  routinely.  This  procedure  is  also  used 
in  case  of  a  time  reversal:  the  last  full  window  before 
the  reversal  and  the  first  full  window  after  the  start 
of  correct  time  progression  are  used  before  the  analysis 
continues.  The  program  indicates  on  the  output  bulletin 
when  data  gaps  or  time  reversals  were  detected,  and 
gives  the  corrective  procedure  used. 


The  fast  f-k  analysis  algorithm 

The  power  at  a  given  frequency  and  wavenumber  is 
computed  as: 


P(f,  k) 


1 

FT 


l  A  (f)exp[2Tri(<f>  (f)-k«r  )] 
n=l  1  11  11 


(1) 
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where 


f  =  frequency 

k  =  vector  wavenumber 

N  =  number  of  data  channels 

n  =  sensor  or  channel  index 

r  =  vector  location  of  the  n'th  sensor  with 

— n 

respect  to  an  arbitrary  origin. 

An(f)exp  [2iri(j)nCf)  1  =  Fourier  transform  of 
the  n'th  data  channel. 

An(f)  is  the  amplitude  part  of  the  transform 
exp[2iri<f>n(f)]  is  the  phase  part  of  the  trans¬ 
form  in  which  <t>n(f)  is  the  phase  angle. 

Expression  (1)  is  evaluated  for  a  matrix  of  wave- 
number  values  at  a  series  of  discrete  frequencies,  as 
specified  in  the  input  parameters;  it  can  be  considered 
as  a  three-dimensional  transform  space  with  frequency 
being  one  dimension  and  the  vector  wavenumber  k  being 
the  other  two  dimensions.  For  computation,  k  is  resolved 
into  a  Cartesian  coordinate  system,  with  k^.  related  to 
geographic  north  and  k  related  to  geographic  east. 

A  wavenumber  value,  say  kg,  is  related  to  the 
phase  velocity  V  by: 

V  -  f/k0 

Thus  phase  velocity  is  inversely  proportional  to  the 

distance  from  the  frequency  axis.  The  locus  of  constant 

values  of  V  is  a  cone  in  f-k  space,  with  the  apex  at 

the  point  f  -  k  =  k  =  0. 

a  y 

For  f-k  analysis  we  calculate  the  power  at  a  matrix 
of  wavenumber  values  separated  by  a  grid  interval  Ak. 
This  is  greatly  speeded  up  by  using  the  following 


(2) 
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relation: 


An(£)exp[2TTiOn(f)  -  (k+Ak)  *  rn)] 

=  An(£)exp[2Tri((j)n(£)  -  k.r^)  ]exp[-2iriA^r  ] 


(3) 


Thus  i£  a  set  of  N  terms  had  been  calculated  for  the 

first  wavenumber  value  kj,  the  values  at  k2  =  k^  +  Ak 

are  obtained  merely  by  multiplying  those  terms  by  a 

factor  exp(-2iri  k«rn) .  Hence  if  a  regular  grid  is  used, 

only  one  set  of  kernels  exp(- 2^1^.  r^)  need  be  generated, 

the  remaining  values  being  obtained  with  successive 

multiplication  by  the  invariant  kernels  exp(-27ri  Ak»r  1 

- nJ  * 

Initially  a  coarse  wavenumber  grid  is  used  to  obtain 
an  approximate  location  for  the  power  maximum;  we  discuss 
the  selection  of  spacing  below.  Once  a  rough  location 
for  the  first  maximum  has  been  obtained,  an  'uphill, 
walk'  approach  is  used  to  determine  the  exact  location 
of  the  peak.  Beginning  at  the  point  of  greatest  power 
on  the  coarse  grid,  the  program  steps  out  in  a  plane 
of  constant  frequency  along  each  of  the  four  cartesian 
coordinate  directions  to  determine  the  direction  in  which 
the  power  is  rising,  and  then  continues  to  compute 
successive  points  in  that  direction  as  long  as  the  power 
is  rising.  When  the  power  begins  to  fall  off  in  the 
direction  being  explored,  a  new  direction  is  determined 
and  followed,  and  the  process  is  repeated  until  a  place 
is  reached  where  the  four  adjacent  points  in  f-k  space 
all  show  lower  power.  The  grid  spacing  is  then  reduced 
Dy  a  factor  of  6  and  the  same  procedure  is  repeated  to 
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refine  the  location  of  the  power  maximum.  The  amount  of 
computation  required  is  about  an  order  of  magnitude  less 
than  would  be  required  for  computing  and  searching  the 
complete  two-dimensional  spectral  matrix. 

All  two-dimensional  peaks  located  in  this  manner 
are  then  checked  to  see  if  they  are  also  maxima  in  the 
frequency  direction  as  well;  such  peaks  are  defined  by: 

3P  _  3P  3P  n 

sr  -  nr  =  tit  =  0 

A  y 

Of  course,  the  extrema  considered  are  just  the  maxima; 
minima  are  not  of  interest. 

The  program  checks  through  the  different  frequencies 
as  follows:  consider  the  power  in  the  wavenumber  plane 
for  the  n'th  frequency.  The  power  maximum  in  this  layer 
is  compared  with  the  power  and  location  of  the  maximum 
in  the  (n-l)'th  and  (n+IJ'th  layers;  for  clarity  of 
exposition  we  call  the  power  at  these  peaks  P  ,  P  , 

and  Pn+1°  If  both  Pn+1  and  Pn-1  are  smaller  than  P^, 
then  Pn  is  flagged  as  a  three-dimensional  maximum.  If 

either  Pn+1  or  P^  is  greater  than  Pn,  a  check  is  made 
on  the  respective  positions  of  the  peaks  in  wavenumber 
space;  suppose  that  Pn_^  >  The  vector  wavenumber 
separation  between  the  location  of  ?n  and  the  location 

of  Pn-1  is  calculated,  and  the  following  logic  is 
followed: 

(a)  The  separation  is  less  than  half  the  width  of 
the  main  lobe  of  the  array  response.  In  this  case,  the 
two  peaks  are  presumed  to  be  part  of  the  same  signal, 
and  P^  is  not  designated  as  a  three-dimensional  maximum. 

(b)  The  separation  is  greater  than  one  halfwidth 
of  the  main  lobe  of  the  array  response.  In  this  case, 


P  if  presumed  to  be  part  of  another  signal.  To  see 

whether  P  itself  is  nevertheless  a  three-dimensional 

maximum,  the  power  is  checked  in  the  (n-l)’th  and  (n+D'th 

layers  at  the  same  vector  wavenumber  as  the  peak  P_n- 

suppose  these  are  Qn_1  and  Qnfl  respectively.  Then  Pn  is 

designated  as  a  three-dimensional  peak  if  PR  >  Qn-1  and 

p  >  o  , .  This  identification  process  is  continued  for 
n  sn+l 

all  the  frequency  layers. 

Wavenumber  peaks  which  are  two-dimensional  bu*  not 
three-dimensional,  i.e.,  those  for  which 


3P  _  9P  _  n 


3P 


i  o 


may  nevertheless  appear  in  the  bulletin,  listed  separately 
from  the  three-dimensional  peaks.  As  a  check  on  the 
validity  of  the  two-dimensional  peaks,  a  group  velocity 
check  is  made,  as  follows.  The  values  of  Af/Ak  are 
calculated  for  each  peak  and  listed  in  the  bulletin, 
where  Af  is  the  frequency  difference  between  adjacent 
layers  and  Ak  is  the  difference  in  wavenumber  position 
of  the  maxima  in  the  two  adjacent  layers.  The  ratio 
Af/Ak  is  a  first-order  approximation  to  the  group 
velocity  of  a  propagating  signal;  thus,  if  two-dimensional 
peak  is  representative  of  the  power  in  a  propagating 
signal  at  that  given  frequency,  the  approximate  group 
velocity  should  have  a  reasonable  value.  A  common  cause 
of  spurious  two-dimensional  peaks  is  leakage  of  power 
(due  to  the  finite  length  of  the  time  window)  from  a 
large  peak  in  the  spectrum.  We  have  shown  previously 
(Smart,  1971)  that  this  leakage  occurs  along  lines  of 
constant  wavenumber  (i.e.,  Ak  !J  0),  so  the  calculated 
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approximation  to  group  velocity  Af/Ak  will  fall  outside 
the  normal  range  of  group  velocities  expected  for  seismic 
signals  (2  to  3  kilometer/second). 

At  each  frequency,  only  the  largest  two-dimensional 
maximum  is  considered  at  first,  in  order  to  avoid  picking 
power  entering  the  array  through  sidclobes  in  the  array 
response .  After  each  peak  is  examined  and  either  accepted 
or  rejected  as  a  signal  detection,  FKCOMB  has  the  option 
of  returning  to  the  frequency  section  being  examined,  in 
order  to  consider  secondary  maxima.  These  smaller  relative 
maxima  may  well  be  sidelobes  of  the  leargest  peak,  and 
in  fact,  the  whole  section  will  be  contaminated  with 
power  from  that  peak.  Thus  it  is  necessary  to  remove  the 
largest  peak,  i.e.,  'strip'  it  out,  before  proceeding 
with  the  anlaysis. 

Thus,  after  a  maximum  in  the  n'th  frequency  section 
has  been  examined,  it  is  removed  by  the  stripping  process 
(see  below),  along  with  its  sidelobes,  and  the  plane  is 
searched  for  another  power  maximum.  This  peak  is  then 
subjected  to  the  previously  described  comparison  with 
the  (n-l)'th  and  (n+l)'th  layers  both  before  and  after 
they  have  themselves  been  stripped.  If  a  detection  is 
made  after  stripping,  the  fact  is  so  indicated  on  the 
output  bulletin.  FKCOMB  at  present  picks  only  the 
primary  and  one  secondary  peak  at  e3ch  frequency. 
Experience  in  using  FKCOMB  to  analyze  long-period  seismic 
data  has  shown  that  signals  are  usually  detected  from  the 
three-dimensional  maxima,  and  the  main  interest  in  the 
two-dimensional  peaks  is  to  provide  a  more  complete 
spectrum  of  the  detected  signal  waveform. 
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The  stripping  process 

The  removal  of  the  first  signal  peak  and  its  associ¬ 
ated  sidclobcs  is  achieved  by  subtracting  an  estimate  of 
an  appropriately  phase-shifted  complex  signal  from  the 
complex  Fourier  transforms  of  the  individual  channels. 

If  a  particular  peak  under  consideration  has  a 
maximum  at  wavenumber  k ' ,  the  estimate  of  the  signal 
at  frequency  f  used  in  the  stripping  process  is: 


Jr  1  A  (f)cxp [2ni (6  (f)  -  k'  *r  )]  (4) 

"  n-1  n  n  n 


This  estimate  is  removed  from  the  Fourier  transform 
set  by  subtracting  it  from  each  of  the  transforms  phase- 
shifted  to  k'  and  then  shifting  the  remnant  back  to  the 
origin  in  the  wavenumber  plane.  The  operation  for  the 
j*th  data  channel  is  thus: 


(AjtOcxpUniUjCf)  -  -  IT  l1AnCf)cxp[2.i(*n(f) 

-  k’ *1^)  J }  cxp(2iri  k*«rj) 


(5) 


From  these  N  filtered  transforms  a  new  wavenumber 
spectrum  is  computed,  and  its  largest  maximum  is  taken 
as  the  second  pick  at  the  frequency  under  consideration, 
to  be  accepted  or  rejected  using  the  same  criteria  as 
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were  applied  to  the  first  pick. 

This  filtering  scheme  is  based  on  two  assumptions: 
first,  that  the  signal  is  a  plane  wave,  although  slight 
departures  from  this  condition  are  not  significant, 
especially  when  the  wavelength  is  of  the  same  order  as 
the  array  diameter.  Second,  that  the  spatial  cnvclopo 
for  the  propagating  waves  is  approximately  flat  over 
distances  equal  to  the  size  of  the  array,  so  that  the 
spatial  spectral  smoothing  function  is  simply  the  array 
response  to  a  plane  wave  whose  phase  velocity  is  infinite 
and  whose  waveform  is  an  impulse. 

The  latter  assumption  bears  on  the  procedure  of 
subtracting  the  mean  Fourier  transform  from  each  of  the 
N  individual  channel  frequency  transforms,  for  in  doing 
that  we  arc  in  effect  subtracting  the  complex  array 
response  from  the  wavenumber  transform,  end  thus  removing 
not  only  the  signal  peak  but  also  all  its  associated 
sidclobcs. 

The  second  wavenumber  spectrum  c>;n  be  expressed  as: 

h  N 

(f3cxp(2iri  (f)  -  k'  *ij)  J 

]  N 

*  FT  J/n^expUHC^Cf)  -  k ••rn)J! 

,2 

•  cxp( 2nik* Je;:p(-2irik«£j  ]  J 


(6) 
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or  alternatively : 


^  l  Aj  tf)cxp[2ni(<>j(0  -  li'Ij)) 

-  J  lll'v  lyncxv[2«H*n(<)  - 

•  cxp( 2nik’  •  r_]  }exp(-2itik^  i^j  j  * 


(7) 


The  squared  modulus  of  the  difference  between  the  two 
complex  terms  between  the  vertical  bars  is  the  estimate 
of  the  wavenumber  spectrum  that  would  be  calculated  if 
the  signal  at  k1  were  not  present. 

The  term  to  the  left  of  the  minus  sign  is  the  original 
unfiltcrcd  complex  wavenumber  spectrum.  The  expression 
in  square  brackets  to  the  right  is  a  constant  transform, 
invariant  from  channel  to  channel.  Thus  the  entire  term 
to  the  right  of  the  minus  sign  is  simply  a  complex  array 
response  multiplied  by  the  amplitude 


|i  l  A  (f) exp l 2ai (♦„(£)  -  L'-I*)!! 

'N  n-1  n 

of  the  signal  estimate  at  f.  .'»cncy  and  phase-shifted 
by  an  amount 
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which  is  the  phase  angle  of  the  signal  estimate,  all 
offset  n  amount  Jc'  from  the  origin  in  the  wavenumber 
plane. 

The  assumption  that  the  spatial  envelope  is  flat 
is  reasonable  when  considering  long-period  surface  waves 
as  recorded  at  ALPA,  LASA,  and  NORSAR,  because  at  the 
periods  of  interest  there  arc  at  most  only  two  cycles  of 
the  signal  within  the  dimension  of  the  array  for  any 
data  window.  Variations  in  instrument  response  will, 
however,  introduce  an  unknown  contribution  to  the  spatial 
spectral  smoothing  function  and  thus  impose  a  limit  on 
the  capability  of  this  stripping  process  to  separato 
signals  of  very  different  amplitude.  This  limit  would, 
of  course,  apply  to  any  array  processing  technique  which 
depended  on  having  constant  or  at  least  known  instrument 
responses . 

An  example  of  the  application  of  the  stripping 
process  is  shown  in  Figure  A.  The  program  FKCONB  does 
not  at  present  furnish  contoured  sections,  but  these 
are  easily  calculated,  and  are  useful  for  illustrating 
the  f-k  analysis. 

The  particular  example  used  here  shows  how  a 
dominant  coherent  peak  is  stripped  away  along  with  its 
sidclobcs,  and  on  recomputing  the  spectrum  in  the  same 


w  ivonumbcr  plane,  the  smaller  second  peak  emerges  clearly. 

S* gnol- to-noisc  ratio  and  F  statistic 

The  signal-to-noisc  ratio  is  defined  as: 


S/N 


P(f,k) 


Av(f )  -  P  (  f ,  k ) 


where  P(f,k)  is  the  power  at  the  maximum  peak,  and  is 
taken  as  the  signal  power  estimate  at  frequency  f.  A v ( f ) 
is  the  average  channel  power  (the  power  averaged  over 
channels)  at  frequency  f: 


Av(f) 


1 

FT 


N 

l 

n»l 


A„(0 


This  value  is  recalculated  after  stripping. 
The  F  statistic  is  defined  as: 


,  P(f’^ 

F  "  Ay(f)  -  P(f.k) 


(N-l) 


(Shumway ,  1971).  An  interesting  aspect  of  the  F  statistic 
is  that  it  allows  us  to  distinguish  energy  arriving  in 
sidclobfes  of  the  array  fro*  energy  in  the  «ain  lobe. 

For  a  *orc  detail  discussion  of  the  application  of  the 
F  statistic  to  seismic  signal  detection  see  Blandford 
(1971);  examples  of  the  application  to  infrasonic 
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array  data  processing  arc  given  by  Smart  and  Fllnn  (1971). 


Selection  of  grid  spacing 

As  mentioned  above,  a  coarrc  wavenumber  grid  is 
chosen  initially  in  order  to  keep  computing  time  to  a 
minimum.  However,  the  spacing  must  be  fine  enough  to 
ensure  that  at  least  one  point  falls  on  a  main  lobe  above 
the  height  of  the  largest  sidclobc  of  the  array  response. 
For  example,  if  the  largest  sidclobc  of  the  array  response 
has  a  peak  9  dB  down  from  the  maximum  of  the  response, 
then  the  grid  spacing  has  to  be  selected  so  that  at  least 
one  point  falls  within  the  9  dB  contour  of  the  main  lobe 
for  any  signal  in  the  velocity  range  of  interest. 

In  the  ease  of  arrays  where  spatial  aliasing  is  a 
problem  (such  as  ALFA)  the  above  condition  must  be 
qualified  by  specifying  the  largest  sidclobc  which  may 
be  expected  to  occur  in  the  velocity  and  frequency  range 
specified  by  the  input  parameters. 

In  the  initial  coarsc-grid  computation  of  the  wave- 
number  spectrum  FKCOMB  employes  an  ir.oscclcs  triangular 
grid  because  it  requires  fewer  points  and  less  computa¬ 
tion  time  than  docs  a  square  grid.  Nevertheless,  the 
convention  has  been  adopted  that  users  will  specify  the 
optimum  square  grid  interval  required  for  the  given 
sensor  array  and  FKCOMB  will  convert  that  to  the 
corresponding  isosceles  triangular  grid  spacing.  Both 
grid  intervals,  DKX  and  DKY,  must  be  specified  and  they 
arc  equal. 
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APPENDIX 


FKCOMB  used  as  the  SAAC  LP  processor 

As  was  mentioned  in  the  introduction  FKCOMB  forms 
the  heart  of  the  SAAC  LP  processing  package.  The  current 
version  can  be  used  directly  on  the  following  kinds  of 
data. 

1.  LASA,  NORSAR  and  ALPA  long  period  data  from  both 
the  SAAC  lowrate  and  highrate  tapes. 

2.  LACA  short  period  from  the  highrate  tapes. 

3.  TFSO  long  period  data. 

The  input  parameters  are  as  follows: 

1)  MOTION,  FOLLOW  -  FORMAT  (2/8) 

Enter  VERTICAL  for  vertical  motion 

LONGITUDINAL  for  horizontal  radial  motion 
TRANSVERSE  for  horizontal  transverse  motion 
(FOLLOW  merely  takes  up  the  extra  letters). 

2)  SWITCH  -  FORMAT  (15) 

Enter  0  for  processing  without  stripping 
1  for  processing  with  stripping 

3)  INC  -  FORMAT  (15) 

Sample  interval  between  start  of  sequential  data 
blocks.  If  this  figure  equals  the  number  of  samples 
used  in  a  data  window  there  is  no  overlap. 

4)  1ST  -  FORMAT  (15) 

Number  of  the  lowest  frequency 


A-l 

fi 


desired  from  the  set  of  discrete  frequencies  (numbers 
increase  with  frequency) . 

5)  ISP  -  FORMAT  (15) 

Number  of  the  highest  frequency  required,, 

NB  The  number  of  any  particular  frequency  is  calculated 
from  the  window  length  and  the  sampling  frate. 

6)  KCOUNT  -  FORMAT  (15) 

Number  of  data  windows  to  be  analyzed. 

7)  SLOP  -  FORMAT  (F10.4) 

Variance  limit  used  in  the  editing  state. 

8)  UNDER,  TOP  -  FORMAT  (2F10.4) 

Minimum  and  maximum  velocities  of  interest.  FKCOMB 
only  searches  the  k-plane  inside  the  lowest  velocity 
limit  and  only  lists  detections  made  between  the  limits. 

9)  GLOWER,  GUPPER  -  FORMAT  (2F10.4) 

Lower  and  upper  group  velocity  limits  defining 
range  over  which  two  dimensional  peaks  are  accepted. 

10)  FSMIN  -  FORMAT  (F10;4) 

Minimum  Fisher  statistic  below  which  no  peaks  are 
listed  as  detections.  Defined  in  the  text. 

11)  DKX.DKY  -  FORMAT  (2F10.6) 

Wavenumber  grid  separation  values,  see  text  for 
discussion  on  how  to  compute  these. 

12)  ANGLE  -  FORMAT  (F10.4) 

Sets  the  inclination  of  the  axis  of  the  cone  to 
be  searched  in  f-k  space.  Normally  ANGLE  is  set  to 
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zero  and  the  cone  axis  is  the  frequency  axis, 

13)  NDAY,  NHR,  NMIN,  NSEC,  NTENTH,  MDAY,  MHR,  MMIN ,  MSEC, 
MTENTH,  ID,  NYEAR  -  FORMAT  (515,  IOX,  515,  3X,  A2,  15) 

NDAY, . . . , NTENTH,  is  the  day,  hour,  minute,  second 
and  decisecond  of  the  start  of  the  first  data  block. 

MDAY, . . . , MTENTH,  defines  the  end  of  the  first  block. 

All  subsequent  blocks  will  have  the  same  length  as  the 
first. 

ID  defines  the  record  type: 

AL:  ALPA 

LA:  LASA  short  period 
LL:  LASA  long  period 
XW:  NORSAR 
TF:  TFO 

NYEAR  is  the  calendar  year  in  which  the  data  were  recorded. 
The  year  is  not  included  in  the  data  records  and  so  this 
parameter  does  not  affect  the  time  search  procedures. 

14)  NSMPLS ,  NCHLS,  LEADER,  NWORDS,  NORMAL,  NDP  -  FORMAT  (6110) 

NSMPLS  =  number  of  sample  times  per  data  record. 

NCHLS  =  number  of  channels  the  user  wishes  to 
retrieve. 

LEADER  =  number  of  2-byte  words  in  the  record  before 
the  first  data  word. 

NWORDS  =  number  of  2-byte  words  which  one  must  pass 
in  going  from  a  given  channel  in  one  frame 
to  the  rest 

NORMAL  =  the  expected  number  of  deciseconds  between 
the  start  of  sequential  data  records. 

NDP  =  5  if  ID  =  XW,  0  otherwise.  NSMPLS,  LEADER, 
NWORDS,  NORMAL  and  NDP  are  invariant  for  a 
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particular  ID.  The  required  values  are  as  follows: 


NSMPLS 

LEADER 

NWORDS 

NORMA 

NDP 

ID 

5 

75 

409 

5 

0 

LA 

10 

20 

131 

100 

0 

LL 

15 

10 

95 

150 

0 

AL 

10 

22 

153 

100 

5 

XW 

60 

6 

36 

600 

0 

TF 

15)  IWORD(I),  ARRAY (I) ,  SENSOR(I) ,  Y(I),  X(I),  I=1,NCHLS 
-  FORMAT  (1X,I3,3X,A4,18X,A2,28X,F7.3,7X,F7.3) 

IWORD(I)  =  sequence  number  of  the  It^1  channel  in 
a  data  frame 

ARRAY(I) ,SENSOR(I) =  are  not  used  by  the  program  but  only 

serve  to  identify  the  array  and  sensor 
to  the  user 

Y (I)  =  north  coordinate  of  the  Ith  sensor 

X(I)  =  east  coordinate  of  the  I  sensor. 
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Figure  2.  Schematic  representation 


SET  S/N  DETECTION 
THRESHOLD.  SET 
VELOCITY  RANGE 


_ 1 _ 

REMOVE  MAXIMUM  PEAKS 
AND  ASSOCIATED  SIDE 
LOBES  IN  F  -  K  SPACE 
PICK  SECOND  MAXIMA 


the  FKCOMB  processor 


222222°°o-o 

ooooooooooo 

Ui|UUj|UU|u(WIWUi|iiiij 

O  O  ^  C  Q  pj  m  4k  — - 

£  £  2  S' ^  N  •'l<n  N  m  « 

ooooddddddo 


*  £  _  ^  H  *4  ft- 

l^ffOfyOOO<N>o 

mcSo«o«Oininmo 

I  •**  «y 


Of»o>o<DrwOOOoc 

omooo«.^o*owJ 

I  I  fU 


sssssssg 


a 

ft 


r w 

9 


a 

oe 

fr¬ 

iz* 


OOOOOOOOOOO 

UlWIMU«UiU(U/UMUUUi 

«ftlf>0ft>Kft)^ft)Of"*<W 

OOOOOOOOOOO 


o  o  © 

U>  ft*  Ui 


*  «r  m 

IZ>  U>  I/* 


a 

a 

o 


o  *  ®  f*  < 


tN  mhOO 


> 

«■ 

o 


z 

u 


o 

z 

C/ 


NtlAWWNO  fr*  ft*  <c  <c 


******  fM  * 


'ft-  O'  «# 

1  — ®  «x 


a 

o 


a 

c- 


ooirrgoirp»l^y^f>. 
*  *  m*'*'*r'r'DD(T\t* 


<oi»(»nro«Nffww 

oir  onw^Mrto^'J 


Jjnjwnnnn 

oooooooo 

Ui  u;  u;  tu  iu  Ui  u<  Ui 

Q  ^  ^  ^  pp  ^  ^ 

OOOOOOOO 


C3 


-^OOOftlft**)* 


ft 

o 


o  «#  aMiMTON 

-Nir  o  O  +  r*  + 
*4  m 


x 

v 

o«r>oo-mcjfw 

4  S  IT  Ohm^o 


O^f'N^KOO 

■*  *«  *  i*  o  o  *2  r 


i 

o 

4/ 

«# 

*» 

o 


*  0  O  0  N  •<  m  • 
^  0  ¥>  r*  m  i 


* 

S 


w  o 

ft  ~ 
S  fr- 


»c  i 


o  o 
o  o 


o  **  o  o 

**<**«*rft» 

(So(Sc({oofS 


ft-  o 

-*  P* 

O  m 


O  O  o  O  o  —  OmOO 

OOOOOOOO WO 

iuuiu/u/uiuiu/iiiu/u, 

5S22=CJ2g; 

dodddddddd 


«2:9°°ooo 

ooooooooo 

u<  IK  lli  IW  Ui  iu  IU  ^ 
T*  *  *  NO 

£  fo^«rsir<(f« 

vrr^or^cv'^o 

**•••*■«■ 

OOOOOOOOO 


So5o7KooCSSiJ 

^OmOHKOOO-^N# 

«*  ft 

*So«ooomoo5** 

c5VoVo#c5r,v0*o*j;r: 


o  a 

ft  a 


a  or 

h-  ►- 

»/*  * 


c 

11/ 

ft 

ft 


ft  ft 
ft  ft 


o 

11/ 

ft 


i 

«• 

x 


ft 

o 


^  ^  ^  W  N  ^  PM  PM  pm  Pj  «  M  ft. 

ooooooooooooo 

JWttiUiWUili'UJIWUiUilWUiUj 

St?*'*0®*'*®*'*? 

OOOOOOOOOOOOO 


O  **  *  *ft  ft  I*. 


'  —  o 


♦  H«f|flWWp<Nft  *  •  * 

fttTIWlWfyfyfyfyfy^^.^® 


o-j  wir  OMrHfv^tr^p, 
WONf  O  *>  r  ft 
•M  P4  ff>  ^  fy  ry  f»,  py 


OOMIOM 


O 

u 


ft 

o 


>irr»«oKK«:irr*fypMa 

JomftJo-oo*^*:.: 


xrwmrMfwrwrwf*  fM 

ooooooooo 

UiUUiU.UiOiUiU.U/ 

~  2  !T  *  ^  ^  «-•  ** 

•  . . 

OOOOOOOOO 


c 

z 

c 


ft-  O  f  *  K  0  <  4  o  M 
•  •  ••••••• 

•v  N  K  h-  ^  y  P  4)K 


o 

U/ 

B 

X  ®  i r  ir  r  *r  4>o 

(*» 


OOKtfNO^rwa*, 

•-0  **  O  ft  ft  «  ,r  «• 

«#•••.  ft  ®  f*  O  ft  ft? 


Ninvc«<r«pM|/ 

•  •  •  • . 

o-^a  OMMifttr  o 


o  o 

P4  M 


:  ft  «  ft- 
•HftH 


®  r»  O 


•■•  5  2  2  5?  *  2  2  2  *  _ _ _ 

V  >  S  ,#Nnnn«  * ^ 

A  «  »  - 

o  N 


* 

M  c/ 

•  ■ — 
1  K 

V> 


wi  I 
*  u. 


£4ipbon//i  «  v 

O  M  1 


T) 

a> 

a 

a 


V) 


"O 

h  < 
o  so 
J  c 


<u 

-c 

H 


Ou 

a 


•  4J 
4-»  V) 
3 

a  ^ 

+j  o 
3  +J 
O  <4-1 
« 
"O 

o  t> 

♦J  T3 
V)  « 
•H  6 
rH 

V) 
4>  W 
i 

•M 

C 

44  O 
O-H 
+-> 
<u  o 

iH  4) 

a+j 

B  a) 

«  TJ 

X 

4)  4> 
-C 
C  *j 
< 

44 
•  *H 

ro 

V) 

4>  U 

C8 

3  4) 

40  Q. 

•h  cl 

U.  CD 


I  1*1  , - “ - WIHIWHMMIIIIII 

*j  ■  IIMIKt  UtlllMi 

Figure  4a.  Frequency  wavenumber  representation  of  time 

vi^?^--ThCkWa'rCnUabfr  planc  15  ^own  the  left  anS 
version  is  shown  on  the  riohr. 
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version  of  the  wavenumber  plane  shown  on  the  left. 
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